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Abstract 

We construct a system of nonequilibrium entropy limiters for the lattice Boltz- 
mann methods (LBM). These limiters erase spurious oscillations without blurring 
of shocks, and do not affect smooth solutions. In general, they do the same work 
for LBM as flux limiters do for finite differences, finite volumes and finite elements 
methods, but for LBM the main idea behind the construction of nonequilibrium en- 
tropy limiter schemes is to transform a field of a scalar quantity — nonequilibrium 
entropy. There are two families of limiters: (i) based on restriction of nonequilibrium 
entropy (entropy "trimming" ) and (ii) based on filtering of nonequilibrium entropy 
(entropy filtering). The physical properties of LBM provide some additional bene- 
fits: the control of entropy production and accurate estimate of introduced artificial 
dissipation are possible. The constructed limiters are tested on classical numeri- 
cal examples: ID athermal shock tubes with an initial density ratio 1:2 and the 
2D lid-driven cavity for Reynolds numbers Re between 2000 and 7500 on a coarse 
100 X 100 grid. All limiter constructions are applicable for both entropic and non- 
entropic quasiequilibria. 
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1 Introduction 



In 1959, S.K. Godunov [T7] demonstrated that a (linear) scheme for a PDE 
could not, at the same time, be monotone and second order accurate. Hence, 
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we should choose between spurious oscillation in high order non-monotone 
schemes and additional dissipation in first order schemes. Flux limiter schemes 
are invented to combine high resolution schemes in areas with smooth fields 
and first order schemes in areas with sharp gradients. 

The idea of fiux limiters can be illustrated by computation of the fiux Fq i 
of the conserved quantity u between a cell marked by and one of two its 
neighbour cells marked by ±1: 

Fo,i = (l-</'(r))/fr + </'W/oT' 

where 'f are low and high resolution scheme fiuxes, respectively, r = 

{uq — — Mo), and 0(r) > is a fiux limiter function. For r close to 1, 

the fiux limiter function 0(r) should be also close to 1. 

Many fiux limiter schemes have been invented during the last two decades [43] . 
No particular limiter works well for all problems, and a choice is usually made 
on a trial and error basis. 



Below are several examples of fiux limiter functions: 
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lod, [36]): 

< 2) (Osher, [10]); 
2)] (monotonised central [42]): 
2)] (superbee, [36]): 
(1 < /3 < 2) (Sweby, [40j). 



The lattice Boltzmann method has been proposed as a discretization of Boltz- 
mann's kinetic equation and is now in wide use in fiuid dynamics and beyond 
(for an introduction and review see [38]). Instead of fields of moments M, the 
lattice Boltzmann method operates with fields of discrete distributions /. This 
allows us to construct very simple limiters that do not depend on slopes or 
gradients. 

All the limiters we construct are based on the representation of distributions 
/ in the form: 

/ = r + ii/-rii^^^, 

where /* is the correspondent quasiequilibrium (conditional equilibrium) for 
given moments M, / — /* is the nonequilibrium "part" of the distribution, 
which is represented in the form "norm x direction" and ||/ — /*|| is the norm 
of that nonequilibrium component (usually this is the entropic norm). Lim- 
iters change the norm of the nonequilibrium component / — /*, but do not 
touch its direction or the equilibrium. In particular, limiters do not change the 
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macroscopic variables, because moments for / and /* coincide. All limiters we 
use are transformations of the form 

/^r + 0x(/-r) (1) 

with > 0. If / — /* is too big, then the limiter should decrease its norm. 

The outline of the paper is as follows. In Sec. [2] we introduce the notions and 
notations from lattice Boltzmann theory we need, in Sec. [3] we elaborate the 
idea of entropic limiters in more detail and construct several nonequilibrium 
entropy limiters for LBM, in Sec. H] some numerical experiments are described: 

(1) ID athermal shock tube examples; 

(2) steady state vortex centre locations and observation of first Hopf bifur- 
cation in 2D lid-driven cavity flow. 

Concluding remarks are given in Sec. [51 



2 Background 

The essence of lattice Boltzmann methods was formulated by S. Succi in the 
following maxim: "Nonlinearity is local, non-locality is linear"[l]. We should 
even strengthen this statement. Non- locality (a) is linear; (b) is exactly and 
explicitly solvable for all time steps; (c) space discretization is an exact oper- 
ation. 

The lattice Boltzmann method is a discrete velocity method. The finite set 
of velocity vectors {vi} {i = l,...m) is selected, and a fluid is described by 
associating, with each velocity Vi, a single-particle distribution function /, = 
fi{x,t) which is evolved by advection and interaction (collision) on a fixed 
computational lattice. The values fi are named populations. If we look at all 
lattice Boltzmann models, one finds that there are two steps: free flight for 
time 6t and a local collision operation. 

The free flight transformation for continuous space is 

fi{x, t + 6t) = fi{x - Vi6t, t). 

After the free flight step the collision step follows: 

Mx) ^ F,i{f,ix)}), (2) 

^ S. Succi, "Lattice Boltzmann at all-scales: from turbulence to DNA transloca- 
tion" , Mathematical Modelling Centre Distinguished Lecture, University of Leices- 
ter, Leicester UK, 15**^ November 2006. 
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or in the vector form 

f{x)^F{f{x)). 

Here, the collision operator F is the set of functions Fi{{fj}) {i = l,...m). 
Each function Fi depends on all fj (j = 1, ...m): new values of the populations 
/, at a point x are known functions of all previous population values at the 
same point. 

The lattice Boltzmann chain "free flight collision free flight collision 
• • • " can be exactly restricted onto any space lattice which is invariant with 
respect to space shifts of the vectors Vi6t {i = 1, ...m). Indeed, free flight trans- 
forms the population values at sites of the lattice into the population values 
at sites of the same lattice. The collision operator ([2]) acts pointwise at each 
lattice site separately. Much effort has been applied to answer the questions: 
"how does the lattice Boltzmann chain approximate the transport equation for 
the moments M?" , and "how does one construct the lattice Boltzmann model 
for a given macroscopic transport phenomenon?" (a review is presented in 
book [Ml). 

In our paper we propose a universal construction of limiters for all possible 
collision operators, and the detailed construction of Fi{{fj}) is not important 
for this purpose. The only part of this construction we use is the local equilibria 
(sometimes these states are named conditional equilibria, quasiequilibria, or 
even simpler, equilibria). 

The lattice Boltzmann models should describe the macroscopic dynamic, i.e., 
the dynamic of macroscopic variables. The macroscopic variables Mi{x) are 
some linear functions of the population values at the same point: M^{x) = 
Y^im£ifi{x), or in the vector form, M{x) = m{f{x)). The macroscopic vari- 
ables are invariants of collisions: 

E^uh = J2^aF,i{f,}) (or m(/) = m(F(/))). 

i i 

The standard example of the macroscopic variables are hydrodynamic fields 
(density-velocity-energy density): {n,nu, E}{x) := Y.A^^'^i^'^'i 1'^} fii^)- -^^^ 
this is not an obligatory choice. If we would like to solve, by LBM methods, 
the Grad equations p2j or some extended thermodynamic equations [25], we 
should extend the list of moments (but, at the same time, we should be ready 
to introduce more discrete velocities for a proper description of these extended 
moment systems). On the other hand, the athermal lattice Boltzmann models 
with a shortened list of macroscopic variables {n, nu] are very popular. 

The quasiequilibrium is the positive fixed point of the collision operator for 
the given macroscopic variables M. We assume that this point exists, is unique 
and depends smoothly on M. For the quasiequilibrium population vector for 
given M we use the notation /]^, or simply /*, if the correspondent value of 
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M is obvious. We use 11* to denote the equilibration projection operation of 
a distribution / into the corresponding quasiequihbrium state: 

n*(/) = 

For some of the colhsion models an entropic description of equilibrium is pos- 
sible: an entropy density function S{f) is defined and the quasiequihbrium 
point fl.]- is the entropy maximiser for given M [26|l39] . 

As a basic example we shall consider the lattice Bhatnagar-Gross-Krook 
(LBGK) model with overrelaxation (see, e.g., p[T2]|23]l28]l38] ) . The LBGK col- 
lision operator is 

F(/):=n*(/) + (2/3-l)(n*(/)-/), (3) 

where /3 G [0, 1]. For (3 = 0, LBGK colhsions do not change /, for /5 = 1/2 
these collisions act as equilibration (this corresponds to the Ehrenfests' coarse 
graining [15] further developed in p^l9ll20] ). for j3 = 1, LBGK collisions act 
as a point reflection with the center at the quasiequihbrium n*(/). 

It is shown [H] that under some stability conditions and after an initial period 
of relaxation, the simplest LBGK collision with overrelaxation [2^11^ provides 
second order accurate approximation for the macroscopic transport equation 
with viscosity proportional to 6t{l — (3)/ (3. 

The entropic LBGK (ELBM) method [5]l20ll26ll39] differs in the definition 
of ([3]) : for (3 = lit should conserve the entropy, and in general has the following 
form: 

F(/):=(l-/5)/ + /3/, (4) 

where / = (1 — a)f + all*{f). The number a = a(/) is chosen so that the 
constant entropy condition is satisfied: S{f) = S{f). For LBGK ([3]), a = 2. Of 
course, for ELBM the entropic definition of quasiequihbrium should be valid. 

In the low-viscosity regime, LBGK suffers from numerical instabilities which 
readily manifest themselves as local blow-ups and spurious oscillations. 

The LBM experiences the same spurious oscillation problems near sharp gra- 
dients as high order schemes do. The physical properties of the LBM schemes 
allows one to construct new types of limiters: the nonequilibrium entropy lim- 
iters. In general, they do the same work for LBM as flux limiters do for finite 
differences, finite volumes and finite elements methods, but for LBM the main 
idea behind the construction of nonequilibrium entropy limiter schemes is to 
limit a scalar quantity — nonequilibrium entropy (and not the vectors or ten- 
sors of spatial derivatives, as it is for flux limiters). These limiters introduce 
some additional dissipation, but all this dissipation could easily be evaluated 
through analysis of nonequilibrium entropy production. 
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Two examples of such limiters have been recently proposed: the positivity 
rule Pl31ll41j and the Ehrenfests' regularisation [7]. The positivity rule just 
provides positivity of distributions: if a collision step produces negative popu- 
lations, then the positivity rule returns them to the boundary of positivity. In 
the Ehrenfests' regularisation, one selects the k sites with highest nonequilib- 
rium entropy (the difference between entropy of the state / and entropy of the 
corresponding quasiequilibrium state /* at a given space point) that exceed a 
given threshold and equilibrates the state in these sites. 

The positivity rule and Ehrenfests' regularisation provide rare, intense and 
localised corrections. It is easy and also computationally cheap to organise 
more gentle transformation with smooth shift of highly nonequilibrium states 
to quasiequilibrium. The following regularisation transformation distributes 
its action smoothly: we can just choose in ([1]) = 0(A5'(/)) with sufficiently 
smooth function (j){AS{f)). Here / is the state at some site, /* is the corre- 
sponding quasiequilibrium state, 5* is entropy, and AS{f) := S{f*) — S{f). 

The next step in the development of the nonequilibrium entropy limiters is in 
the usage of local entropy filters. The filter of choice here is the median filter: it 
does not erase sharp fronts, and is much more robust than convolution filters. 

An important problem is: "how does one create nonequilibrium entropy lim- 
iters for LBM with non-entropic quasiequilibria?" . We propose a solution 
of this problem based on the nonequilibrium KuUback entropy. For entropic 
quasiequilibrium the Kullback entropy approach gives the same entropic lim- 
iters. In thermodynamics, Kullback entropy belongs to the family of Massieu- 
Planck-Kramers functions (canonical or grandcanonical potentials). 



3 Nonequilibrium entropy limiters for LBM 

3.1 Positivity rule 

There is a simple recipe for positivity preservation Pl31|41j : to substitute 
nonpositive lQ{f){x) by the closest nonnegative state that belongs to the 
straight line 

|A/(x) + (l-A)n*(/(x))| A gm} (5) 

defined by the two points, f{x) and corresponding quasiequilibrium. This op- 
eration is to be applied pointwise, at points of the lattice where positivity 
is violated. The coefficient A depends on x too. Let us call this recipe the 
positivity rule (Fig. [1]). This recipe preserves positivity of populations and 
probabilities, but can affect accuracy of approximation. The same rule is nec- 
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Fig. 1. Positivity rule in action. The motions stops at the positivity boundary. 

essary for ELBM (jlj) when the positive "mirror state" / with the same entropy 
as / does not exists on the straight hne (jS]). 

3.2 Ehrenfests' regularisation 

To discuss methods with additional dissipation, the entropic approach is very 
convenient. Let entropy S{f) be defined for each population vector / = (/«) 
(below we use the same letter S for local in space entropy, and hope that 
context will make this notation always clear). We assume that the global 
entropy is a sum of local entropies for all sites. The local nonequilibrium 
entropy is 

AS{f):=S{n-S{f), (6) 
where /* is the corresponding local quasiequilibrium at the same point. 

The Ehrenfests' regularisation [6|7] provides "entropy trimming": we moni- 
tor local deviation of / from the corresponding quasiequilibrium, and when 
AS{f){x) exceeds a pre-specified threshold value 6, perform local Ehrenfests' 
steps to the corresponding quasiequilibrium: / i— /* at those points. 

So that the Ehrenfests' steps are not allowed to degrade the accuracy of LBGK 
it is pertinent to select the k sites with highest AS > 6. The a posteriori 
estimates of added dissipation could easily be performed by analysis of entropy 
production in Ehrenfests' steps. Numerical experiments show (see, e.g., [6,7J) 
that even a small number of such steps drastically improve stability. 

To avoid the change of accuracy order "on average" , the number of sites with 
this step should be < 0{Nh/L) where N is the total number of sites, h is 
the step of the space discretization and L is the macroscopic characteristic 
length. But this rough estimate of accuracy in average might be destroyed 
by concentration of Ehrenfests' steps in the most nonequilibrium areas, for 
example, in the boundary layer. In that case, instead of the total number of 
sites N in 0{Nh/L) we should take the number of sites in a specific region. 
The effects of concentration could be easily analysed a posteriori. 
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3.3 Smooth limiters of nonequilibrium entropy 



The positivity rule and Ehrenfests' regularisation provide rare, intense and 
localised corrections. Of course, it is easy and also computationally cheap to 
organise more gentle transformation with a smooth shift of highly nonequilib- 
rium states to quasiequilibrium. The following regularisation transformation 
distributes its action smoothly: 

/^r + </.(A5(/))(/-r). (7) 

The choice of function (p is highly ambiguous, for example, = 1/(1 + aAS'^) 
for some a > and A; > 0. There are two significantly different choices: (i) 
ensemble- independent (p (i.e., the value of depends on local value of AS* 
only) and (ii) ensemble-dependent 0, for example 

1 + (Ag/(aE(Ag)))^-V^ 
"^^^^^ - 1 + iAS/iaEiASW ' 

where E(AS') is the average value of AS" in the computational area, k > 1, 
and a > 1. For small A^, (f){AS) ^ 1 and for A^ > aE(A^), (f){AS) tends 



to y q;E(AS')/AS'. It is easy to select an ensemble-dependent with control 
of total additional dissipation. 



3.4 Monitoring of total dissipation 



For given (3, the entropy production in one LBGK step in quadratic approxi- 
mation for AS* is: 

5lbgk^ ~ [1 - (2/3 - If ] ^ A5(cc), 

X 

where x is the grid point, AS{x) is nonequilibrium entropy at point x, 
^LBGK'S' is the total entropy production in a single LBGK step. It would be 
desirable if the total entropy production for the limiter Sii^S was small relative 

to ^LBGK'S': 

5iin,S < SqSlbgkS. (9) 

A simple ensemble-dependent limiter (perhaps, the simplest one) for a given 
5o operates as follows. Let us collect the histogram of the AS{x) distribution, 
and estimate the distribution density, p{AS). We have to estimate a value 
ASq that satisfies the following equation: 

/ p{AS){AS-ASo)dAS = 5o[l-{2p-lY] p{AS)ASdAS. (10) 

JASo Jo 
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In order not to affect distributions with small expectation of AS, we choose 
a threshold ASt = max{ASo,6}, where 6 is some predefined value (as in 
the Ehrenfests' regularisation). For states at sites with AS" > ASt "we pro- 



vide homothety with quasiequihbrium center /* and coefficient y ASt/ AS (in 
quadratic approximation for nonequilibrium entropy): 



3. 5 Median entropy filter 

The limiters described above provide pointwise correction of nonequilibrium 
entropy at the "most nonequilibrium" points. Due to the pointwise nature, 
the technique does not introduce any nonisotropic effects, and provides some 
other benefits. But if we involve the local structure, we can correct local non- 
monotone irregularities without touching regular fragments. For example, we 
can discuss monotone increase or decrease of nonequilibrium entropy as regular 
fragments and concentrate our efforts on reduction of "speckle noise" or "salt 
and pepper noise". This approach allows us to use the accessible resource of 
entropy change ([9]) more thriftily. 

Among all possible filters, we suggest the median filter. The median is a more 
robust average than the mean (or the weighted mean) and so a single very 
unrepresentative value in a neighborhood will not affect the median value 
significantly. Hence, we suppose that the median entropy filter will work better 
than entropy convolution filters. 

The median filter considers each site in turn and looks at its nearby neighbours. 
It replaces the nonequilibrium entropy value AS at the point with the median 
of those values AS'med, then updates / by the transformation (ITT]) with the 



homothety coefficient y A5'mod/A5'. The median, A^med, is calculated by first 
sorting all the values from the surrounding neighbourhood into numerical order 
and then replacing that being considered with the middle value. For example, 
if a point has 3 nearest neighbors including itself, then after sorting we have 
3 values AS*: ASi < AS2 < AS^,- The median value is A^med = AS2- For 9 
nearest neighbors (including itself) we have after sorting AS'med = AS^. For 
27 nearest neighbors AS'med = AS'14. 

We accept only dissipative corrections (those resulting in a decrease of AS', 
AS'med < AS) because of the second law of thermodynamics. The analogue 
of (ITO!) is also useful for acceptance of the most significant corrections. 

Median filtering is a common step in image processing [31] for the smoothing 
of signals and the suppression of impulse noise with preservation of edges. 




(11) 
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3.6 Entropic steps for non-entropic quasiequilibria 

Beyond the quadratic approximation for nonequilibrium entropy all the logic of 
the above mentioned constructions remain the same. There exists only one sig- 
nificant change: instead of a simple homothety flTT]) with coefficient ^ ASt/ AS 
the transformation ([7]) should be applied, where the multiplier is a solution 
of the nonlinear equation 

sir+(j)if-n) = s{n-As,. 

This is essentially the same equation that appears in the definition of ELBM 
steps dl). 

More differences emerge for LBM with non-entropic quasiequilibria. The main 
idea here is to reason that non-entropic quasiequilibria appear only because of 
technical reasons, and approximate continuous physical entropic quasiequilib- 
ria. This is not an approximation of a density function, but an approximation 
of measure, i.e., from the cubature formula: 

i 
i 

The discrete populations fi are connected to continuous (and sufficiently 
smooth) densities f{v) by cubature weights fi ~ Wif{vi). These weights for 
quasiequilibria are found by moment and flux matching conditions |37j. It 
is impossible to approximate the BGS entropy / /In/dt? just by discretiza- 
tion (to change integration by summation, and continuous distribution / by 
discrete /i), because cubature weights appear as additional variables. Never- 
theless, the approximate discretization of the KuUback entropy S'k [3U] does 
not change its form: 

^k(/) = - / f{v) In dt; ^ - E /an (-|) , (12) 

because fi/ f* approximates the ratio of functions f{v)/f*{v) and J2ifi--- 
gives the integral / f{v) . . .dv approximation. Here, in f|T2l) . the state /* is the 
quasiequilibrium with the same values of the macroscopic variables as /. More- 
over, for given values of the macroscopic variables, S^if) achieves its maxi- 
mum at the point f = f* (both for continuous and for discrete distributions). 
The corresponding maximal value is zero. Below, S'k is the discrete Kullback 
entropy. If the approximate discrete quasiequilibrium /* is non-entropic, we 
can use — S'k(/) instead of AS{f). 

For entropic quasiequilibria with perfect entropy the discrete Kullback entropy 
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gives the same AS": — 5'k(/) = AS{f). Let the discrete entropy have the 
standard form for an ideal (perfect) mixture [27]. 



After the classical work of Zeldovich [H], this function is recognised as a 
useful instrument for the analysis of kinetic equations (especially in chemical 
kinetics [21]). If we define /* as the conditional entropy maximum for given 
Mj = Y.k^jkfk, then 

j 

where fij{M) are the Lagrange multipliers (or "potentials"). For this entropy 
and conditional equilibrium we find 

A^ = ^(r)-^(/) = ^/an(^A^, (13) 

if / and /* have the same moments, m(/) = m(/*). The right hand side 
of m is -5k(/). 

In thermodynamics, the Kullback entropy belongs to the family of Massieu- 
Planck-Kramers functions (canonical or grandcanonical potentials). There is 
another sense of this quantity: Sk is the relative entropy of / with respect to 
/* [T81I35] . 

In quadratic approximation, 

-^K(/) = E/an(|)-E^^^^- 



3. 7 ELBM collisions as a smooth limiter 



On the base of numerical tests, the authors of [41] claim that the positivity 
rule provides the same results (in the sense of stability and absence/presence 
of spurious oscillations) as the ELBM models, but ELBM provides better 
accuracy. 

For the formal definition of ELBM (jl]) our tests do not support claims that 
ELBM erases spurious oscillations (see below). Similar observation for Burgers 
equation was previously published in [1]. We understand this situation in the 
following way. The entropic method consists at least of three components: 

(1) entropic quasiequilibrium, defined by entropy maximisation; 
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(2) entropy balanced collisions (jl]) that have to provide proper entropy bal- 
ance; 

(3) a method for the solution of the transcendental equation S{f) = S{f) to 
find a = a{f) in 

It appears that the first two items do not affect spurious oscillations at all, 
if we solve the equation for a{f) with high accuracy. Additional viscosity 
is, potentially, added by explicit analytic formulas for a{f). In order not to 
decrease entropy, errors in these formulas always increase dissipation. This 
can be interpreted as a hidden transformation of the form ([7]), where the 
coefficients in depend also on /*. 

3.8 Monotonic and double monotonic limiters 

Two monotonicity properties are important in the theory of nonequilibrium 
entropy limiters: 

(1) a limiter should move the distribution to equilibrium: in all cases of ([T]) 
< < 1. This is the dissipativity condition which means that limiters 
never produce negative entropy. 

(2) a limiter should not change the order of states on the line: if for two 
distributions with the same moments, / and /', AS{f) > AS{f') before 
the limiter transformation, then the same inequality should hold after the 
limiter transformation too. For example, for the limiter ([7]) it means that 
AS{f* + x(f){AS{f* + x{f - f*)){f - /*)) is a monotonically increasing 
function of x > 0. 

In quadratic approximation, 

A5(r + x(/-r)) = a;2A5(/), 
AS{r + x^Sir + x{f - /*))(/ - D) = x'<p\x'AS{f)), 

and the second monotonicity condition transforms into the following require- 
ment: y(j){y'^s) is a monotonically increasing (not decreasing) function of y > 
for any s > 0. 

If a limiter satisfies both monotonicity conditions, we call it "double mono- 
tonic" . For example, Ehrenfests' regularisation satisfies the first monotonicity 
condition, but obviously violates the second one. The limiter ([H]) violates the 
first condition for small AS, but is dissipative and satisfies the second one in 
quadratic approximation for large AS". The limiter with = 1/(1 + aAS'') al- 
ways satisfies the first monotonicity condition, violates the second if > 1/2, 
and is double monotonic (in quadratic approximation for the second condi- 
tion), ifO<A;<l/2. The threshold limiters (II ip are also double monotonic. 
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Of course, it is not forbidden to use any type of limiters under the local and 
global control of dissipation, but double monotonic limiters provide some nat- 
ural properties automatically, without additional care. 



4 Numerical experiment 

To conclude this paper we report some numerical experiments conducted to 
demonstrate the performance of some of the proposed nonequilibrium entropy 
limiters for LBM from Sec. [31 



4.1 Velocities and quasiequilibria 



We will perform simulations using both entropic and non-entropic quasiequi- 
libria, but we always work with an athermal LBM model. Whenever we use 
non-entropic quasiequihbria we employ Kullback entropy f|T3|) . 

In ID, we use a lattice with spacing and time step 6t = 1 and a discrete 
velocity set {vi, V2, V3} := {0, —1, 1} so that the model consists of static, left- 
and right-moving populations only. The subscript i denotes population (not 
lattice site number) and /i, /2 and /s denote the static, left- and right-moving 
populations, respectively. The entropy is 5* = —H, with 



(see, e.g., [27]) and, for this entropy, the local entropic quasiequilibrium state 
/* is available explicitly: 



H = h log(/i/4) + /2 log(/2) + h log(/3) 



fl 




(14) 



P 
6 



((3n+ 1) -2V1 + 3m2), 



where 



1 



(15) 



P 



u := — 



P 
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The standard non-entropic polynomial quasiequilibria [38] are: 



/* = ^(l-3n + 3n2) 
o 

/* = ^(1 + 3^ + 3^2) 





(16) 



In 2D, we employ a uniform 9-speed square lattice with discrete velocities 
{vi\i = 0, 1, . . . 8}: vo = 0, Vi = (cos((i - l)n/2), sm{{i - 1)tt/2)) for i = 
1,2,3,4, Vi = V2(cos((z - 5)f + f),sin((i - 5)f + f)) for i = 5,6,7,8. The 
numbering /o, fi,...,fs are for the static, east, north, west, south, north- 
east, northwest, southwest and southeast-moving populations, respectively. 
As usual, the entropic quasiequilibrium state, /*, can be uniquely determined 
by maximising an entropy functional 

5(/) = -E/aog(^), 

subject to the constraints of conservation of mass and momentum [2]: 

Here, the lattice weights, Wi, are given lattice-specific constants: Wo = 4/9, 
W^i,2,3,4 = 1/9 and PV5,6,7,8 = 1/36. Analogously to ( fT5l) . the macroscopic vari- 
ables p and u = (mi,M2) are the zeroth and first moments of the distribution 
/, respectively. The standard non-entropic polynomial quasiequilibria [38] are: 

/• = pw-.(i + 3... + H(i|^i)!_!|!). (18) 

4.2 LBGK and ELBM 

The governing equations for LBGK are 

U{x + v,,t + l) = f:{x,t) + {2f3 - l)(/;(x, t) -Mx,t)), (19) 
where (3 = l/{2u + l). 
For ELBM (jlj) the governing equations are: 

fi{x + t;„ t + 1) = (1 - (3)f:{x, t) + (3Mx, t), (20) 



14 



with /? as above and / = {l — a)f + af*. The parameter, a, is chosen to satisfy 
a constant entropy condition. This involves finding the nontrivial root of the 
equation 

Siil-a)f + an = S{f). (21) 

To solve (I2T!) numerically we employ a robust routine based on bisection. The 
root is solved to an accuracy of 10~^^ and we always ensure that the returned 
value of a does not lead to a numerical entropy decrease. We stipulate that 
if, at some site, no nontrivial root of ( |2T1) exists we will employ the positivity 
rule instead (Fig. [1]). 

4-3 Shock tube 

The ID shock tube for a compressible athermal fiuid is a standard benchmark 
test for hydrodynamic codes. Our computational domain will be the interval 
[0, 1] and we discretize this interval with 801 uniformly spaced lattice sites. 
We choose the initial density ratio as 1:2 so that for x < 400 we set p = 1.0 
else we set p = 0.5. We will fix the kinematic viscosity of the fiuid at u = 10~^. 

4.3.1 Comparison of LBGK and ELBM 

In Fig. [2] we compare the shock tube density profile obtained with LBGK 
(using entropic quasiequilibria f|T^ ) and ELBM. On the same panel we also 
display both the total entropy S{t) := J2xS{x,t) and total nonequilibrium 
entropy AS{t) := J2x ^S{x,t) time histories. As expected, by construction, 
we observe that total entropy is (effectively) constant for ELBM. On the other 
hand, LBGK behaves non-entropically for this problem. In both cases we ob- 
serve that nonequilibrium entropy grows with time. 

As we can see, the choice between the two collision formulas LBGK ( IT9!) 
or ELBM (l20l) does not affect spurious oscillation, and reported regularisa- 
tion [29j is, perhaps, the result of approximate analytical solution of the equa- 
tion (!?!]) . Inaccuracy in the solution of can be interpreted as a hidden 
nonequilibrium entropy limiter. But it should be mentioned that the entropic 
method consists not only of the collision formula, but, what is important, in- 
cludes a special choice of quasiequilibrium that could improve stability (see, 
e.g., [13]). Indeed, when we compare ELBM with LBGK using either eutopic or 
standard polynomial quasiequilibria, there appears to be some gain in employ- 
ing entropic quasiequilibria (Fig. [3]). We observe that the post-shock region 
for the LBGK simulations is more oscillatory when polynomial quasiequilibria 
are used. In Fig. [3] we have also included a panel with the simulation result- 
ing from a much higher viscosity {u = 3.3333 x 10~^). Here, we observe no 
appreciable differences in the results of LBGK and ELBM. 
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Fig. 2. Density and profile of the 1:2 athermal shock tube simulation with = 10 
after 400 time steps using (a) LBGK ([19]); (b) ELBM ([20]). In this example, no 
negative population are produced by any of the methods so the positivity rule is 
redundant. For ELBM in this example, (j2ip always has a nontrivial root. Total 
entropy and nonequilibrium entropy time histories are shown in panels (c), (d) and 
(e), (f) for LBGK and ELBM, respectively. 




Fig. 3. Density and velocity profile of the 1:2 isothermal shock tube simula- 
tion after 400 time steps using (a) LBGK (jl9p with polynomial quasiequilib- 
ria (HI]) [i/ = 3.3333 x 10"^]. (j^) LBGK ([H]) with entropic quasiequilibria <^ 
[v = 3.3333 X 10-2]; (c) ELBM ^ [v = 3.3333 x lO'^]; (d) LBGK ([H]) with 
polynomial quasiequilibria ()16p [u = 10^^]; (e) LBGK (jlOp with entropic quasiequi- 
libria ^[u = 10-9]; (f) ELBM ^ [v = lO^^]. 
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4-3.2 Nonequilibrium entropy limiters. 



Now, we would like to demonstrate just a representative sample of the many 
possibilities of limiters suggested in Sec. [31 In each case the limiter is im- 
plemented by a post-processing routine immediately following the collision 
step (either LBGK or ELBM ([20])). Here, we will only consider LBGK 
collisions and entropic quasiequilibria f[T^ . 

The post-processing step adjusts / by the update formula: 

f^r + ^{As){f-n, 

where A^* is defined by ([6]) and is a limiter function. 
For the Ehrenfests' regularisation one would choose 



<P{AS)ix) 



AS{x) < 6, 
otherwise, 



where 5 is a pre-specified threshold value. Furthermore, it is pertinent to select 
just k sites with highest AS > S. This limiter has been previously applied to 
the shock tube problem in [6][7[[8] and we will not reproduce those results here. 

Instead, our first example will be the following smooth limiter: 

For this limiter, we will fix /c = 1/2 (so that the limiter is double monotonic in 
quadratic approximation to entropy) and compare the density profiles for a = 
6/ (E(A5')^), 6 = 0.1, 0.01, 0.001. We have also ensured an ensemble-dependent 
limiter because of the dependence of a on the average E(A5'). As with Fig. [21 
we accompany each panel with the total entropy and nonequilibrium entropy 
histories. Note the different scales for nonequilibrium entropy. Note also that 
entropy (necessarily) now grows due to the additional dissipation. 

Our next example (Fig. [5]) considers the threshold filter ([TOl) . In this example 
we choose the estimates A^o = 5E(AS), lOE(A^), 20E(A5) and fix the tol- 
erance (5 = so that the influence of the threshold alone can be studied. Only 
entropic adjustments are accepted in the limiter: ASt < AS. As the threshold 
increases, nonequilibrium entropy grows faster and spurious begin to appear. 

Finally, we test the median filter (Fig. [6]). We choose a minimal filter so that 
only the nearest neighbours are considered. As with the threshold filter, we 
introduce a tolerance 6 and we try the values S = 10~^, 10""^, 10~^. Only 
entropic adjustments are accepted in the limiter: AS'med < AS. 
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Fig. 4. Density and profile of the 1:2 athermal shock tube simulation with = 10 ^ 
after 400 time steps using LBGK (fTU|) and the smooth limiter with k = 1/2, 
a = 5/(E(AS')'=) and (a) 5 = 0.1; (b) 5 = 0.01 and (c) 5 = 0.001. Total entropy and 
nonequilibrium entropy time histories for each parameter set {k, a{5)} are displayed 
in the adjacent panels. 

We have seen that each of the examples we have considered (Fig. HI Fig. [5] 
and Fig. [6]) is capable of subduing spurious post-shock oscillations compared 
with LBGK (or ELBM) on this problem (cf. Fig. [2]). Of course, by limiting 
nonequilibrium entropy the result is necessarily an increase in entropy. 

From our experiences our recommendation is that the median filter is the 
superior choice amongst all the limiters suggested in Sec. [31 The action of the 
median filter is found to be both extremely gentle and, at the same time, very 
effective. 



4.4 Lid- driven cavity 



Our second numerical example is the classical 2D lid-driven cavity flow. A 
square cavity of side length L is filled with fluid with kinematic viscosity v 
(initially at rest) and driven by the cavity lid moving at a constant velocity 
[uq, 0) (from left to right in our geometry). 
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Fig. 5. Density and profile of the 1:2 athermal shock tube simulation with = 10 ^ 
after 400 time steps using LBGK (|19|) and the threshold limiter (|10p with (a) 
i^St = 5E(A5); (b) ASt = 10E(A5) and (c) A^* = 20E(A5). Total entropy and 
nonequilibrium entropy time histories for each threshold AS't are displayed in the 
adjacent panels. 

We will simulate the flow on a 100 x 100 grid using LBGK regularised with 
the median filter limiter. Unless otherwise stated, we use entropic quasiequilib- 
ria (1171) . The implementation of the filter is as follows: the filter is not applied 
to boundary nodes; for nodes which immediately neighbour the boundary the 
stencil consists of the 3 nearest neighbours (including itself) closest to the 
boundary; for all other nodes the minimal stencil of 9 nearest neighbours is 
used. 



We have purposefully selected such a coarse grid simulation because it is read- 
ily found that, on this problem, unregularised LGBK fails (blows-up) for all 
but the most modest Reynolds numbers Re := Luq/u. 



4-4-1 Steady-state vortex centres 

For modest Reynolds number the system settles to a steady state in which the 
dominant features are a primary central rotating vortex, with several counter- 
rotating secondary vortices located in the bottom-left, bottom- right (and pos- 
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Fig. 6. Density and profile of the 1:2 athermal shock tube simulation with v = 10~^ 
after 400 time steps using LBGK (I19p and the minimal median limiter with (a) 
5 = 10~^; (b) 5 = 10~^ and (c) 5 = 10"^. Total entropy and nonequilibrium 
entropy time histories for each tolerance 5 are displayed in the adjacent panels. 



sibly top-left) corners. 



Steady state has been extensively investigated in the literature. The study 
of Hou et al [21] simulates the flow over a range of Reynolds numbers using 
unregularised LBGK on a 256 x 256 grid. Primary and secondary vortex centre 
data is provided. We compare this same statistic for the present median filtered 
coarse grid simulation. We will employ the same convergence criteria used 
in ^24j. Namely, we deem that steady state has been reached by ensuring 
that the difference between the maximum value of the stream function for 
successive 10,000 time steps is less that 10^^. The stream function, which is 
not a primary variable in the LBM simulation, is obtained from the velocity 
data by integration using Simpson's rule. Vortex centres are characterised as 
local extrema of the stream function. 

We compare our results with the LBGK simulations in [21] and [H]. To align 
ourselves with these studies we specify the following boundary condition: lid 
profile is constant; remaining cavity walls are subject to the "bounce-back" 
condition [38j. In our simulations, the initial uniform fluid density profile is 
p = 2.7 and the velocity of the lid is = 1/10 (in lattice units). 
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Collected in Table [H for Re = 2000, 5000 and 7500, are the coordinates of 
the primary and secondary vortex centres using (a) unregularised LBGK; (b) 
LBGK with median filter limiter {6 = 10~^); (c) LBGK with median filter lim- 
iter {S = 10~^), all with non-entropic polynomial quasiequilibria (1181) . Lines 
(d), (e) and (f) are the same but with entropic quasiequilibria f|T7|) . The re- 
maining lines of Table [1] are as follows: (g) literature data [23] (unregularised 
LBGK on a 256 x 256 grid); (h) literature data [?T] (positivity rule); (i) litera- 
ture data [41j (ELBM). With the exception of (g), all simulation are conducted 
on a 100 X 100 grid. The top-left vortex does not appear at Re = 2000 and 
no data was provided for it in [41J at Re = 5000. The unregularised LBGK 
Re = 7500 simulation blows-up in finite time and the simulation becomes 
meaningless. The y-coordinate of the two lower- vortices at Re = 5000 in (i) 
appear anomalously small and were not reproduced by our experiments with 
the positivity rule (not shown). 

We have conducted two runs of the experiment with the median filter param- 
eter 6 = 10~^ and 6 = 10""^. Despite the increased number of realisations the 
vortex centre locations remain effectively unchanged and we detect no signif- 
icant variation between the two runs. This demonstrates the gentle nature of 
the median filter. At Reynolds Re = 2000 the median filter has no effect at all 
on the vortex centres compared with LBGK. 

We find no significant differences between the experiments with entropic and 
non-entropic polynomial quasiequilibria in this test. 

The coordinates of the primary vortex centre for unregularised LBGK at Re = 
5000 are already quite inaccurate as LBGK begins to lose stability. Stability 
is lost entirely at some critical Reynolds number 5000 < Re < 7500 and the 
simulation blows-up. 

Furthermore, we have agreement (within grid resolution) with the data given 
in [23]. Also compiled in Table [1] is the data from the limiter experiments 
conducted in [41j (although not explicitly discussed in the language of limiters 
by the authors of that work). In [41j the authors give vortex centre data for 
the positivity rule (Fig. [T]) and for ELBM (which we interpret as containing a 
hidden limiter). In [IT] the positivity rule is called FIX-UP. 

As Reynolds number increases the flow in the cavity is no longer steady and a 
more complicated flow pattern emerges. On the way to a fully developed tur- 
bulent flow, the lid-driven cavity flow is known to undergo a series of period 
doubling Hopf bifurcations. On our coarse grid, we observe that the coordi- 
nates of the primary vortex centre (maximum of the stream function) is a very 
robust feature of the flow, with little change between coordinates (no change 
in y-coordinates) computed at Re = 5000 and Re = 7500 with the median fil- 
ter. On one hand, because of this observation it becomes inconclusive whether 
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Table 1 

Primary and secondary vortex centre coordinates for the lid-driven cavity flow at 
Re = 2000, 5000, 7500. 

Primary Lower-left Lower-right Top-left 



Rc .(• jj x ij :r ij :r ij 



2000 


(a) 


0.5253 


0.5455 


0.0909 


0.1010 


0.8384 


0.1010 


Not applicable 


2000 


(b) 


0.5253 


0.5455 


0.0909 


0.1010 


0.8384 


0.1010 


Not applicable 


2000 


(c) 


0.5253 


0.5455 


0.0909 


0.1010 


0.8384 


0.1010 


Not applicable 


2000 


(d) 


0.5253 


0.5455 


0.0909 


0.1010 


0.8384 


0.1010 


Not applicable 


2000 


(e) 


0.5253 


0.5455 


0.0909 


0.1010 


0.8384 


0.1010 


Not applicable 


2000 


(f) 


0.5253 


0.5455 


0.0909 


0.1010 


0.8384 


0.1010 


Not applicable 


2000 


(g) 


0.5255 


0.5490 


0.0902 


0.1059 


0.8471 


0.0980 


Not applicable 


2000 


(h) 


0.5200 


0.5450 


0.0900 


0.1000 


0.8300 


0.0950 


Not applicable 


2000 


(i) 


0.5200 


0.5500 


0.0890 


0.1000 


0.8300 


0.1000 


Not applicable 


5000 


(a) 


0.5152 


0.6061 


0.0808 


0.1313 


0.7980 


0.0707 


0.0505 0.8990 


5000 


(b) 


0.5152 


0.5354 


0.0808 


0.1313 


0.8081 


0.0808 


0.0606 0.8990 


5000 


(c) 


0.5152 


0.5354 


0.0808 


0.1313 


0.8081 


0.0808 


0.0707 0.8889 


5000 


(d) 


0.5152 


0.5960 


0.0808 


0.1313 


0.8081 


0.0808 


0.0505 0.8990 


5000 


(e) 


0.5152 


0.5354 


0.0808 


0.1313 


0.8081 


0.0808 


0.0606 0.8990 


5000 


(f) 


0.5152 


0.5354 


0.0808 


0.1313 


0.8081 


0.0808 


0.0707 0.8889 


5000 


(g) 


0.5176 


0.5373 


0.0784 


0.1373 


0.8078 


0.0745 


0.0667 0.9059 


5000 


(h) 


0.5150 


0.5680 


0.0950 


0.0100 


0.8450 


0.0100 


Not available 


5000 


(i) 


0.5150 


0.5400 


0.0780 


0.1350 


0.8050 


0.0750 


Not available 


7500 
7500 


(a) 
(b) 


0.5051 


0.5354 


0.0707 


0.1515 


0.7879 


0.0707 


0.0606 0.8990 


7500 


(c) 


0.5051 


0.5354 


0.0707 


0.1515 


0.7879 


0.0707 


0.0707 0.8889 


7500 


(d) 
















7500 


(e) 


0.5051 


0.5354 


0.0707 


0.1515 


0.7879 


0.0707 


0.0606 0.8990 


7500 


(f) 


0.5051 


0.5354 


0.0707 


0.1515 


0.7879 


0.0707 


0.0707 0.8889 


7500 


(g) 


0.5176 


0.5333 


0.0706 


0.1529 


0.7922 


0.0667 


0.0706 0.9098 
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the median limiter is adding too much additional dissipation. On the other 
hand, a more studious choice of control criteria may indicate that the first 
bifurcation has already occurred by Re = 7500. 

4-4-2 First Hopf bifurcation 

A survey of available literature reveals that the precise value of Re at which 
the first Hopf bifurcation occurs is somewhat contentious, with most current 
studies (all of which are for incompressible flow) ranging from around Re = 
7400-8500 Pl32ll33j . Here, we do not intend to give a precise value because 
it is a well observed grid effect that the critical Reynolds number increases 
(shifts to the right) with refinement (see, e.g.. Fig. 3 in [33] )• Rather, we 
will be content to localise the first bifurcation and, in doing so, demonstrate 
that limiters are capable of regularising without effecting fundamental flow 
features. 

To localise the first bifurcation we take the following algorithmic approach. 
Entropic quasiequilibria are in use. The initial uniform fluid density profile 
is p = 1.0 and the velocity of the lid is Uq = 1/10 (in lattice units). We 
record the unsteady velocity data at a single control point with coordinates 
(L/16, 13L/16) and run the simulation for 5000 non-dimensionless time units 
(5000L/mo time steps). Let us denote the final 1% of this signal by (Msig,fsig)- 
We then compute the energy (£2-norm normalised by non-dimensional 
signal duration) of the deviation of Ugig from its mean: 



(23) 

where |tisig| and denote the length and mean of u^ig, respectively. We 
choose this robust statistic instead of attempting to measure signal amplitude 
because of numerical noise in the LBM simulation. The source of noise in LBM 
is attributed to the existence of an inherently unavoidable neutral stability 
direction in the numerical scheme (see, e.g., [8]). 

We opt not to employ the "bounce-back" boundary condition used in the pre- 
vious steady state study. Instead we will use the diffusive Maxwell boundary 
condition (see, e.g., [HI), which was first applied to LBM in [1]. The essence 
of the condition is that populations reaching a boundary are reflected, propor- 
tional to equilibrium, such that mass-balance (in the bulk) and detail-balance 
are achieved. The boundary condition coincides with "bounce-back" in each 
corner of the cavity. 

To illustrate, immediately following the advection of populations consider the 
situation of a wall, aligned with the lattice, moving with velocity M„aii and 
with outward pointing normal to the wall in the negative ^/-direction (this is 



Uo\U, 



Slg I 
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the situation on the hd of the cavity with M„aii = uq). The implementation 
of the diffusive Maxwell boundary condition at a boundary site (x, y) on this 
wall consists of the update 

+ 1) = 7/;(uwaii), z = 4, 7,8, 

with 

^ f2ix, y, t) + /5(x, y, t) + /6(x, t) 

/4*(Mwall) + /7*(Mwall) + /s ("wall) 

Observe that, because density is a linear factor of the quasiequilibria f|T7|) . 
the density of the wall is inconsequential in the boundary condition and can 
therefore be taken as unity for convenience. As is usual, only those populations 
pointing in to the fluid at a boundary site are updated. Boundary sites do not 
undergo the collisional step that the bulk of the sites are subjected to. 

We prefer the diffusive boundary condition over the often preferred "bounce- 
back" boundary condition with constant lid profile. This is because we have 
experienced difficulty in separating the aforementioned numerical noise from 
the genuine signal at a single control point using "bounce-back" . We remark 
that the diffusive boundary condition does not prevent unregularised LBGK 
from failing at some critical Reynolds number Re > 5000. 

Now, we conduct an experiment and record ( 123|) over a range of Reynolds 
numbers. In each case the median filter limiter is employed with parameter 
5 = 10~^. Since the transition between steady and periodic flow in the lid- 
driven cavity is known to belong to the class of standard Hopf bifurcations 
we are assured that oc Re |16j. Fitting a line of best fit to the resulting 
data localises the first bifurcation in the lid-driven cavity flow to Re = 7135 
(Fig. [7]). This value is within the tolerance of Re = 7402 ± 4% given in [33] for 
a 100 X 100 grid. We also provide a (time averaged) phase space trajectory and 
Fourier spectrum for Re = 7375 at the monitoring point (Fig. [8] and Fig. [9]) 
which clearly indicate that the first bifurcation has been observed. 



5 Conclusions 

Entropy and thermodynamics are important for stability of the lattice Boltz- 
mann methods. It is now clear: after almost 10 years of work since the pub- 
lication of proved this statement (the main reviews are [5)1281139] ). The 
question is now: "how does one utilise, optimally, entropy and thermody- 
namic structures in lattice Boltzmann methods?" . In our paper we attempt to 
propose a solution (temporary, at least). Our approach is applicable to both 
entropic as well as for non-entropic polynomial quasiequilibria. 
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Fie. 7. Plot of energy squared, El as a function of Reynolds number, Re, using 
LBGK regularised with the median filter limiter with 5 = 10~^ on a 100 x 100 grid. 
Straight lines are lines of best fit. The intersection of the sloping line with the x-axis 
occurs close to Re = 7135. 

We have constructed a system of nonequilibrium entropy limiters for the lattice 
Boltzmann methods (LBM): 

• the positivity rule that provides positivity of distribution; 

• the pointwise entropy limiters based on selection and correction of most 
nonequilibrium values; 

• filters of nonequilibrium entropy, and the median filter as a filter of choice. 

All these limiters exploit physical properties of LBM and allow control of total 
additional entropy production. In general, they do the same work for LBM as 
flux limiters do for finite differences, finite volumes and finite elements meth- 
ods, and come into operation when sharp gradients are present. For smoothly 
changing waves, the limiters do not operate and the spatial derivatives can be 
represented by higher order approximations without introducing non-physical 
oscillations. But there are some differences too: for LBM the main idea behind 
the construction of nonequilibrium entropy limiter schemes is to limit a scalar 
quantity — the nonequilibrium entropy — or to delete the "salt and pepper" 
noise from the field of this quantity. We do not touch the vectors or tensors 
of spatial derivatives, as it is for flux limiters. 

Standard test examples demonstrate that the developed limiters erase spurious 
oscillations without blurring of shocks, and do not affect smooth solutions. The 
limiters we have tested do not produce a noticeable additional dissipation and 
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Fig. 8. Velocity components as a function of time for the signal (?Xsig,fsig) at the 
monitoring point (L/16, 13L/16) using LBGK regularised with the median filter 
limiter with 6 = 10^^ on a 100 x 100 grid (Re = 7375). Dots represent simulation 
results and the solid line is a 100 step time average of the signal. 

allow us to reproduce the first Hopf bifurcation for 2D lid-driven cavity on a 
coarse 100 x 100 grid. At the same time the simplest median filter deletes the 
spurious post-shock oscillations for low viscosity. 

Perhaps, it is impossible to find one best nonequilibrium entropy limiter for 
all problems. It is a special task to construct the optimal limiters for a specific 
classes of problems. 
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Fig. 9. Relative amplitude spectrum for the signal ttgig at the monitoring point 
(I//16, 13L/16) using LBGK regularised with the median filter limiter with 5 = 10~^ 
on a fOO X fOO grid (Re = 7375). We measure a dominant frequency of w = 0.525. 
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